%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This script converts trace element sources, and reformats some data for
% plotting. This second bit can be cleaned up and removed...it was written
% very early, prior to many later changes

%The only thing you need to run this is goodTabNames
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


normalizationLabel='PUM';



% Read in trace element sources and select the one for normalization
[xlsNumbers, xlsText,xlsRAW] = xlsread('AAA_PETROGEN_SourcesAndDs', 'TraceElementOrder');
[ElementRow,ElementColumn]= find(strcmp(xlsRAW,'Elements')==1); 
[LastElementRow,LastElementColumn]= find(strcmp(xlsRAW,'LastElement')==1); 
% Matrix with element labels
targetStrings_Trace = xlsRAW(ElementRow+1:LastElementRow-1, ElementColumn)'; 


if any(strcmp(targetStrings_Trace,PetrogenTraceElement_Strings)<1) ==1
    disp('Opps, trace element order seems to have changed! Please Check. Exiting')
    return
end


% Import matrix with constant partition coefficients. Note we do change U
% and Th based on pressure later in the Petrogen code.
[DRow,DColumn]= find(strcmp(xlsRAW,'FirstTESource')==1);
[LastDRow,LastDColumn]= find(strcmp(xlsRAW,'LastTESource')==1);

AllTESources = xlsRAW(ElementRow+1:LastElementRow-1, DColumn:LastDColumn);
AllTESources = cell2mat(AllTESources);
AllTESources=AllTESources';
AllTESources_Names = xlsRAW(ElementRow, DColumn:LastDColumn);





EuCol = find(strcmp(PetrogenTraceElement_Strings,'Eu'));
SmCol = find(strcmp(PetrogenTraceElement_Strings,'Sm'));
GdCol = find(strcmp(PetrogenTraceElement_Strings,'Gd'));
clear DataRatios
for gg = 1:size(ElementRatios_Strings,1)
    AllTESources_DataRatios(:,gg) = AllTESources(:,ElementRatiosIndex(gg,1)) ./AllTESources(:,ElementRatiosIndex(gg,2));
    AllTESources_DataRatios(:,21) = 2.*(AllTESources(:,EuCol)./0.154)./((AllTESources(:,SmCol)./0.406)+(AllTESources(:,GdCol)./0.544));    
end



[a,TraceElements4Normalization_index]=ismember(normalizationLabel,AllTESources_Names); 
TraceElements4Normalization = AllTESources(TraceElements4Normalization_index,:); 
TraceElements4Normalization_DataRatios = AllTESources_DataRatios(TraceElements4Normalization_index,:); 



[a,TraceElements4Normalization_index]=ismember('WHDMM',AllTESources_Names); 
TraceDMM = AllTESources(TraceElements4Normalization_index,:); 
DataRatios_TraceDMM = AllTESources_DataRatios(TraceElements4Normalization_index,:); 


[a,TraceElements4Normalization_index]=ismember('PUM',AllTESources_Names); 
TracePUM = AllTESources(TraceElements4Normalization_index,:); 
DataRatios_TracePUM= AllTESources_DataRatios(TraceElements4Normalization_index,:); 


for kk = 1:size(ElementRatios_Strings,1)
    RatioLabels{kk}=sprintf('%s/%s',ElementRatios_Strings{kk,1}, ElementRatios_Strings{kk,2});
end

% PetrogenTraceElement_Strings
% Petrogen_Co
% Petrogen_Co_Name
% Petrogen_BulkCompositionsNormalized
% Petrogen_Dmatrix


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%



% play(Audio);
% DMM_source_decision = input('Do you want to change all default TE sources to something else? \n [Y/N]:','s');
% if strcmp(DMM_source_decision,'Y')
%     DMM_source_decision = input('Change all models to: [PUM/WHDMM/TSDMM]?:','s');
%     layerplot_TESource(:)={DMM_source_decision};
%     SourceComp = sprintf('%s TE source',DMM_source);% Set Trace Element Source Composition for labels
% else 
%     SourceComp = sprintf('%s TE source','Default');
% end
% 




 for ifmm =  1:size(goodTabNames,2)

     OldSource_Ratios = eval(sprintf('Petrogen_DataRatios_TESource_%s', goodTabNames{ifmm}));
     OldSource_Abundances = eval(sprintf('Petrogen_TESource_%s', goodTabNames{ifmm}));
     OldSource_Name =eval(sprintf('Petrogen_TESourceName_%s', goodTabNames{ifmm})); 
 
     NewSourceName = layerplot_TESource{ifmm};

     
     [a,TraceElements4Normalization_index]=ismember(NewSourceName,AllTESources_Names);
     if TraceElements4Normalization_index>0
         NewSource_Abundances = AllTESources(TraceElements4Normalization_index,:);
         NewSource_Ratios = AllTESources_DataRatios(TraceElements4Normalization_index,:);
     else
         %if no source is given, default to the one run in Petrogen
         NewSource_Abundances = OldSource_Abundances;
         NewSource_Ratios = OldSource_Ratios;
         layerplot_TESource{ifmm} = OldSource_Name; 
         
     end

     
     DataRatios = bsxfun(@rdivide, ...
         eval(sprintf('DataRatios_%s', goodTabNames{ifmm})),...
         OldSource_Ratios./NewSource_Ratios);
     assignin('base',(sprintf('DataRatios_%s', goodTabNames{ifmm})),DataRatios);
     
     
     DataTrace = bsxfun(@rdivide, ...
         eval(sprintf('%s_Trace', goodTabNames{ifmm})),...
         OldSource_Abundances./NewSource_Abundances);
     assignin('base',(sprintf('%s_Trace', goodTabNames{ifmm})),DataTrace);
     
 end

 names4legend_all_sources = strcat(names4legend_all,'-',layerplot_TESource'); 
 

%%

[emorbs1,emorbs2]=find(Galeetal2013_MajorTrace_DataRatios(:,6)./DataRatios_TracePUM(6)>1.5);
[dmorbs1,dmorbs2]=find(Galeetal2013_MajorTrace_DataRatios(:,6)./DataRatios_TracePUM(6)<0.8);
[Hirsch1,Hirsch2]=find(Galeetal2013_MajorTrace_DataRatios(:,11)./DataRatios_TracePUM(11)>1.03 & Galeetal2013_MajorTrace_DataRatios(:,11)./DataRatios_TracePUM(11)<1.28);
[nmorbs1,nmorbs2]=find(Galeetal2013_MajorTrace_DataRatios(:,6)./DataRatios_TracePUM(6)>=0.8 & Galeetal2013_MajorTrace_DataRatios(:,6)./DataRatios_TracePUM(6)<=1.5);
[PossibleGar1,PossibleGar2]=find(Galeetal2013_MajorTrace_DataRatios(:,11)./DataRatios_TracePUM(11)>3);
%[emorbs1_re,emorbs2_re]=find(Gale_5_Master_TE_Iso_Edited_DataRatios(:,6)./DataRatios_PUM(6)>1.5);



%% collects information for plotting

for pp = 1:size(goodTabNames,2)
    
    
    TabNameW_here = goodTabNames{pp};
    %WGP  = 100.*eval(sprintf('%s_Info(:,13)',TabNameW_here));
    [a,xxxvalue]=ismember('percentGar',ConstraintOptions);
    WGP = 100.*eval(sprintf('%s_InfoALL(:,xxxvalue)',TabNameW_here));
    
    [a,xxxvalue]=ismember('percentGar_onaxis',ConstraintOptions);
    WGP_onaxis = 100.*eval(sprintf('%s_InfoALL(:,xxxvalue)',TabNameW_here));
    
    [a,xxxvalue]=ismember('StabilityField ',ConstraintOptions);
    
    if any((strfind(TabNameW_here,'IM')))==1
        zeta_N  = eval(sprintf('%s_InfoALL(:,%d)',TabNameW_here,xxxvalue));
        zeta_N(zeta_N==1)=1;
        zeta_N(zeta_N~=1)=0;
        WGP = 100.*zeta_N;
        WGP_onaxis=100.*zeta_N;
    end
    
    numBulks = size(eval(TabNameW_here),1)/11;
    Garnet =[1:3];
    Garnet_ALL = [];
    NoGarnet =[4:9];
    NoGarnet_ALL = [];
    for llll = 1:numBulks
        NoGarnet_ALL = [NoGarnet_ALL  (NoGarnet+(llll-1).*11)];
        Garnet_ALL = [Garnet_ALL  (Garnet+(llll-1).*11)];
    end
    
    
    [a,xxxvalue]=ismember('ZetaGar',ConstraintOptions);
    if any((strfind(TabNameW_here,'SM')))==1
        WGP  = 100.*eval(sprintf('%s_InfoALL(:,%d)',TabNameW_here,xxxvalue));
        WGP_onaxis  = 100.*eval(sprintf('%s_InfoALL(:,%d)',TabNameW_here,xxxvalue));

    end
    
    %     if any((strfind(TabNameW_here,'SM')))==1
    %         WGP  = 100.*eval(sprintf('%s_Info(:,13)',TabNameW_here));
    %         WGP(1:3)=[100];
    %         WGP(4:9)=[0];
    %     end
    
    assignin('base',sprintf('%s_WGP',TabNameW_here),WGP);
    assignin('base',sprintf('%s_WGP_onaxis',TabNameW_here),WGP_onaxis);
    
    Crust_FP  = eval(sprintf('%s_Info(:,16)',TabNameW_here));
    assignin('base',sprintf('%s_CrustFP',TabNameW_here),Crust_FP);
    
    
    % [ElementRow,ElementColumn]= find(strcmp(xlsRAW,'total F')==1);
    % [DRow,DColumn]= find(~cellfun(@isempty,regexp(xlsRAW_strings,'PlotOn')) ==1);
    % FN = cell2mat(xlsRAW(ElementRow+1:DRow(end), ElementColumn));
    FN  = 100.*eval(sprintf('%s_Info(:,19)',TabNameW_here));
    assignin('base',sprintf('%s_FN',TabNameW_here),FN);
    
    % [ElementRow,ElementColumn]= find(strcmp(xlsRAW,'F_P')==1);
    % [DRow,DColumn]= find(~cellfun(@isempty,regexp(xlsRAW_strings,'PlotOn')) ==1);
    % FP_Full = cell2mat(xlsRAW(ElementRow+1:DRow(end), ElementColumn));
    FP_Full  = 100.*eval(sprintf('%s_Info(:,17)',TabNameW_here));
    assignin('base',sprintf('%s_FP_Full',TabNameW_here),FP_Full);
    
    U  = eval(sprintf('%s_Info(:,5)',TabNameW_here));
    assignin('base',sprintf('%s_U',TabNameW_here),U);
    
    
    
    %ThExcess  = eval(sprintf('%s_Info(:,5)',TabNameW_here));
    %assignin('base',sprintf('%s_U',TabNameW_here),U);
    
    [a,xxxvalue]=ismember('Water',ConstraintOptions);
    WATER = eval(sprintf('%s_InfoALL(:,xxxvalue)',TabNameW_here));
    assignin('base',sprintf('%s_Water',TabNameW_here),WATER);
    
    
    [a,xxxvalue]=ismember('Tdiff',ConstraintOptions);
    temp = eval(sprintf('%s_InfoALL(:,xxxvalue)',TabNameW_here));
    assignin('base',sprintf('%s_Tdiff',TabNameW_here),temp);
    
    
    
    [a,xxxvalue]=ismember('Gar2SpTransition',ConstraintOptions);
    temp = eval(sprintf('%s_InfoALL(:,xxxvalue)',TabNameW_here));
    assignin('base',sprintf('%s_GSBoundary',TabNameW_here),temp);
    
    
    [a,xxxvalue]=ismember('MeltExtractedPb',ConstraintOptions);
    temp = eval(sprintf('%s_InfoALL(:,xxxvalue)',TabNameW_here));
    assignin('base',sprintf('%s_MeltExtractP',TabNameW_here),temp);
    
    
    [a,xxxvalue]=ismember('LastMeltPb',ConstraintOptions);
    temp = eval(sprintf('%s_InfoALL(:,xxxvalue)',TabNameW_here));
    assignin('base',sprintf('%s_LastMeltP',TabNameW_here),temp);
    
    
    
    
end
